package com.asa.utils;

import java.math.BigDecimal;
import java.math.RoundingMode;
import java.text.NumberFormat;

import com.asa.HanShu;
import com.asa.hutils.MathHUtils;
import com.asa.numanaly.F;

/**
 * 基础数据计算方法
 * 
 * @author Administrator
 * 
 */
public class MathUtils {

	/**
	 * 选择牛顿法是因为直接获取原函数即可生成结果，而且拟合过程还算比较快
	 * 
	 * @param x0
	 *            拟合的初始值
	 * @param hanShu
	 *            提交函数接口，调用原函数
	 * @return 求出零点的结果
	 */
	public static double niudun(double x0, HanShu hanShu) {
		int m = 50;
		double q = 0.000001;
		double e = 0.000001;
		double v = hanShu.hanshu(x0);
		if (Math.abs(v) < e) {
			return x0;
		}
		for (int i = 0; i < m; i++) {
			double x1 = x0 - v / dao(x0, hanShu);
			v = hanShu.hanshu(x1);
			// //System.out.println(i+"  "+x0+"  "+v);
			if (Math.abs(x1 - x0) < q || Math.abs(v) < e) {
				break;
			}
			x0 = x1;
		}
		
		
		
		return x0;
	}
	
	/**
	 * 求函数的导数值
	 * @param x0	你希望的求导点
	 * @param hanShu	你提交的初始函数
	 * @return	返回导数求值
	 */
	public static double dao(double x0, HanShu hanShu) {
		// TODO Auto-generated method stub
		return (hanShu.hanshu(x0 + 0.000001) - hanShu.hanshu(x0)) / 0.000001;
	}
	

	/**
	 * 多元一次方程求解，采用的解决方案是高斯法，一般的高斯法矩阵a的对角线上不能为0（因为被除数不能为0），所以需要以以一个最小值来拟合0
	 * 这样精细化程度有限，但实现是最灵活的
	 * 
	 * @param a n元n次矩阵
	 * @param b	n元矩阵
	 * @return 解集
	 */
	public static double[] duoyuanyici(double a[][],double b[]){
		a = MathHUtils.copy(a);
		for (int i = 0; i < a.length; i++) {
			if (a[i][i]==0) {
				a[i][i] = 0.0000001;
			}
		}
	
		double x[] = new double[b.length];
		double sum = 0;
		double xmult = 0;
		for (int k = 0; k < b.length-1; k++) {
			for (int i = k+1; i < b.length; i++) {
				xmult = a[i][k]/a[k][k];
				a[i][k]=xmult;
				for (int j = k+1; j < b.length; j++) {
					a[i][j]=a[i][j]-(xmult)*a[k][j];
				}
				b[i]=b[i]-(xmult)*b[k];
			}
		}
		x[b.length-1] = b[b.length-1]/a[b.length-1][b.length-1];
		for (int i = b.length-2; i >=0; i--) {
			sum = b[i];
			for (int j = i+1; j < x.length; j++) {
				sum = sum - a[i][j]*x[j];
			}
			x[i] = sum/a[i][i];
		}
		
		
		return x;
	}
	
	

	/**
	 * 多元一次方程组的解，使用的是LR分解的方式求解
	 * @param a
	 * @param b
	 * @return
	 */
	
	public static double[] duoyuanyici_LR(double[][] a,double[] b) {
		
		
//		double[][] a = {{2,1,1},
//				{1,3,2,},
//				{1,2,2}};
				
				
				double[][][] doolittle_LR = F.Doolittle_LR(a);
				
//				double[] b = {4,6,5};
				double[][] L = doolittle_LR[0];
				double[][] R = doolittle_LR[1];
				
				double[] y = new double[a.length];
				double[] x = new double[a.length];
				
				for (int i = 0; i < y.length; i++) {
					double asa = 0;
					for (int k = 0; k < i; k++) {
						asa += L[i][k] * y[k];
					}
					
					y[i] = b[i] - asa;
				}
				
				
				for (int i = x.length-1; i >= 0; i--) {
					double asa = 0;
					for (int k = i+1; k < x.length; k++) {
						asa += R[i][k] * x[k];
					}
					
					x[i] = (y[i] - asa) / R[i][i];;
				}
				
				
//				MathHUtils.print(x);

		
		return x;
		
		
		
		
		
		
		
		
		
		
	}
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	/**
	 * 根据传入的点集，拟合一个通过所有点的函数
	 * @param x	x点集，需要与y一一对应
	 * @param y	y点集，需要与x一一对应
	 * @return 返回拟合出来的方程 y = c0 + c1*x + c2*c^2 + c3*c^3 ...
	 */
	public static double[] hangshuhuoqu(double x[],double y[]){
		int n = y.length;
		double[] c = new double[n];
		c[0] = y [0];
		for (int k = 1; k < c.length; k++) {
			double d = x[k]-x[k-1];
			double u = c[k-1];
			for (int i = k-2; i >= 0; i--) {
				u = u*(x[k]-x[i])+c[i];
				d = d*(x[k]-x[i]);
			}
			c[k] = (y[k]-u)/d;
		}
		return c;
	}
	
	/**
	 * 将函数转化为一个我认识的函数
	 * @param c
	 * @return
	 */
	public static HanShu hangshuhuoqu2hanshu(final double[] c){
		HanShu result = new HanShu() {
			
			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				return 0;
			}
			
			@Override
			public double hanshu(double x) {
				// TODO Auto-generated method stub
				double asa = 0;
				for (int i = 0; i < c.length; i++) {
					asa += c[i] * Math.pow(x, i);
				}
				
				
				return asa;
			}
		};
		
		return result;
	}

	
	
	
	
	
	
	
	
	
	
	/**
	 * 龙贝格算法求函数积分
	 * @param a	积分的初始位置
	 * @param b	积分的结束位置
	 * @param hanShu	你期望求解的函数
	 * @return	积分值
	 */
	public static double longBeiGe(double a,double b, HanShu hanShu) {
		
		int M =20;
		double h = b-a;
		double [][] R = new double[M][M];
		
		for (int n = 1; n < R.length; n++) {
			h = h/2;
			double all=0;
			for (int i = 0; i < Math.pow(2, n-1); i++) {
				all= all+hanShu.hanshu(a+(2*i-1)*h);
			}
			R[n][0] = R[n-1][0]/2 + h*all;
			for (int m = 1; m < n; m++) {
				R[n][m] = R[n][m-1]+(R[n][m-1]-R[n-1][m-1])/(Math.pow(4, m)-1);
			}

		}
//		for (int i = 0; i < R.length; i++) {
//			for (int j = 0; j < R.length; j++) {
//				System.out.print(R[i][j]+"\t");
//			}
//			//System.out.println();
//		}
		return R[R.length-1][R.length-2];
	}
	
	
	/**
	 * 四阶龙格库塔方法求解常微分方程的解
	 * @param t	t的初始值
	 * @param x	对应t的初始值得x值
	 * @param tt	要求的t值对应x
	 * @param hanShu	偏微分方程值
	 * @return	要求的点的x的值
	 */
	public static double siJieLongGe(double t,double x,double tt,HanShu hanShu){
		int m = 256;
//		double t = 1;
//		double x = 2;//x的初始值，即t=1时x=2
		double h = (tt - t)/m;//1+h（步长）*m（次数）就是t的结尾了
		
		//double e = Math.abs(u(t)-x);
		//System.out.println(0+" "+t+" "+x+" ");
		for (int k = 0; k < m; k++) {
			double f1 = h*hanShu.hanshu(t,x);
			double f2 = h*hanShu.hanshu(t+h/2,x+f1/2);
			double f3 = h*hanShu.hanshu(t+h/2,x+f2/2);
			double f4 = h*hanShu.hanshu(t+h,x+f3);
			x = x+(f1+2*f2+2*f3+f4)/6;
			t = t+h;
			//e = Math.abs(u(t)-x);
			////System.out.println(k+" "+t+" "+x+" ");//+e);
		}
		return x;
	}
	
	
	
	
	
	

	
	
	
	
	
	
	
	
	
	/**
	 * 小数点后精度整合
	 * @param a
	 * @param dianhou
	 * @return
	 */
	public static Double jindu(Double a,int dianhou) {
		
		
		BigDecimal setScale = new BigDecimal(a).setScale(dianhou, BigDecimal.ROUND_HALF_UP);

		
		return setScale.doubleValue();
		
	}
	
	
	
	

}
